--------------------------------------------------------------------------------
                       # Installing & loading packages
--------------------------------------------------------------------------------
packages <- c("readxl", "dplyr","tidyr","nlme", "mgcv", "itsadug", "flexOR", 
              "parameters",  "rmcorr","ggplot2",  "grid", "tiff", "gridExtra", "patchwork")
# install.packages(setdiff(packages, rownames(installed.packages())))
lapply(packages, library, character.only = TRUE)

# update.packages(checkBuilt=TRUE, ask=FALSE) 
# installed.packages()[,c('Package','Version','LibPath')]
# remove.packages("rmcorr")

--------------------------------------------------------------------------------
                        # Main results of all models
--------------------------------------------------------------------------------
data <- read_excel("data.xlsx")
names <- c("Subject", "Group", "Music", "Preference") # Music=Tonality
data[,names] <- lapply(data[,names], factor)
dat = data
rm(data)

## when Calculating β & 95% CI, Standardize the continuous variables.
# dat$Rating <- scale(dat$Rating)
# dat$PerceivedValence <- scale(dat$PerceivedValence)
# dat$PerceivedArousal <- scale(dat$PerceivedArousal)
# dat$FeltValence <- scale(dat$FeltValence)
# dat$FeltArousal <- scale(dat$FeltArousal)

dat$SubjByMusic = interaction(dat$Subject, dat$Music)
minst = tapply(dat$StandardizedTime, dat$SubjByMusic, min)
dat$MinStTime = minst[as.character(dat$SubjByMusic)]
dat$AR.start = FALSE
dat[dat$MinStTime == dat$StandardizedTime, "AR.start"]=TRUE
--------------------------------------------------------------------------------
# Base model
M0 = bam(Rating ~ Group + s(StandardizedTime, by=Group) + 
           s(StandardizedTime, Subject, bs="fs", m=1), 
           method="ML", data=dat)
# summary(M0)$p.table
# summary(M0)$s.table
# save(M0, file = "Result_M0.RData")
# M0 <- get(load("Result_M0.RData"))
# summary(M0)

# Get autocorrelations
acf_values <- acf(resid(M0), plot=FALSE)
rho <- acf_values$acf[2]  # Lag 1 value
cat("Estimated rho (AR1 coefficient):", round(rho, 2), "\n")

# Plot model diagnostics
tiff(filename="FigureS1.tiff", width=6, height=4.5, units="in", res=900)
acf(resid(M0), plot=TRUE, xlab="Lag", ylab="ACF", las=1, cex.lab=1, font.lab=2, main=" ")
dev.off()
--------------------------------------------------------------------------------
# Group main effect
M1a = bam(Rating ~ Group + s(StandardizedTime, by=Group) + 
            s(StandardizedTime, Subject, bs="fs", m=1),
            AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M1a)
aicc_M1a <- AICc(M1a)
model_parameters(M1a)

# Plot model diagnostics
check_resid(M1a, split_pred='AR.start', select=3, ask=FALSE)
check_resid(M1a, select=1, ask=FALSE)
k.check(M1a)

# Music main effect
M1b = bam(Rating ~ Music + s(StandardizedTime, by=Music) + 
            s(StandardizedTime, Subject, bs="fs", m=1),
            AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M1b)
aicc_M1b <- AICc(M1b)
model_parameters(M1b)

# Group-Music interaction with smooth interaction 
M1c = bam(Rating ~ Group*Music + s(StandardizedTime, by=interaction(Group,Music)) + 
            s(StandardizedTime, Subject, bs="fs", m=1),
            AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M1c)
aicc_M1c <- AICc(M1c)
model_parameters(M1c)
# Note: M1c=M1 Music=Tonality

# Group-Music interaction without smooth interaction
M1d = bam(Rating ~ Group*Music + s(StandardizedTime) + 
            s(StandardizedTime, Subject, bs="fs", m=1),
            AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M1d)
aicc_M1d <- AICc(M1d)
model_parameters(M1d)
--------------------------------------------------------------------------------
 # Repeated Measures Correlation between perceived and felt Valence and arousal 
--------------------------------------------------------------------------------  
# Calculate rmcorr 
stimuli <- read_excel("C:/Users/Administrator/Desktop/Stimuli.xlsx")  
# stimuli <- read_excel("Stimuli.xlsx")
names <- c("Subject","Group","Music")
stimuli[,names] <- lapply(stimuli[,names], factor)
str(stimuli)
summary(stimuli)

rmc1<-rmcorr(participant= Subject, dataset= stimuli, 
             measure1= PerceivedValence, measure2= PerceivedArousal,
             CI.level=0.95, CIs="bootstrap", nreps=10000, bstrap.out=FALSE) 
rmc1

rmc2<-rmcorr(participant= Subject, dataset= stimuli, 
             measure1= PerceivedValence, measure2= FeltValence,
             CI.level=0.95, CIs="bootstrap", nreps=10000, bstrap.out=FALSE) 
rmc2

rmc3<-rmcorr(participant= Subject, dataset= stimuli, 
             measure1= PerceivedValence, measure2= FeltArousal,
             CI.level=0.95, CIs="bootstrap", nreps=10000, bstrap.out=FALSE) 
rmc3 

rmc4<-rmcorr(participant= Subject, dataset= stimuli, 
             measure1= PerceivedArousal, measure2= FeltValence,
             CI.level=0.95, CIs="bootstrap", nreps=10000, bstrap.out=FALSE) 
rmc4

rmc5<-rmcorr(participant= Subject, dataset= stimuli,
             measure1= PerceivedArousal, measure2= FeltArousal,
             CI.level=0.95, CIs="bootstrap", nreps=10000, bstrap.out=FALSE) 
rmc5

rmc6<-rmcorr(participant= Subject, dataset= stimuli,
             measure1= FeltValence, measure2= FeltArousal,
             CI.level=0.95, CIs="bootstrap", nreps=10000, bstrap.out=FALSE) 
rmc6  
--------------------------------------------------------------------------------
M2a = bam(Rating ~ Group*Music + Preference + PerceivedValence + PerceivedArousal +
            s(StandardizedTime, by=interaction(Group,Music)) + 
            s(StandardizedTime, Subject, bs="fs", m=1),
            AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M2a)
aicc_M2a <- AICc(M2a)
model_parameters(M2a)


M2b = bam(Rating ~ Group*Music + Preference + FeltValence + FeltArousal +
            s(StandardizedTime, by=interaction(Group,Music)) + 
            s(StandardizedTime, Subject, bs="fs", m=1),
            AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M2b)
aicc_M2b <- AICc(M2b)
model_parameters(M2b)

--------------------------------------------------------------------------------
M3a = bam(Rating ~ Group*Music + Preference + PerceivedValence +
            s(StandardizedTime, by=interaction(Group,Music)) + 
            s(StandardizedTime, Subject, bs="fs", m=1),
            AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M3a)
aicc_M3a <- AICc(M3a)
model_parameters(M3a)


M3b = bam(Rating ~ Group*Music + Preference + PerceivedArousal +
            s(StandardizedTime, by=interaction(Group,Music)) + 
            s(StandardizedTime, Subject, bs="fs", m=1),
            AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M3b)
aicc_M3b <- AICc(M3b)
model_parameters(M3b)


M3c = bam(Rating ~ Group*Music + Preference + FeltValence +
            s(StandardizedTime, by=interaction(Group,Music)) + 
            s(StandardizedTime, Subject, bs="fs", m=1),
            AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M3c)
aicc_M3c <- AICc(M3c)
model_parameters(M3c)


M3d = bam(Rating ~ Group*Music + Preference + FeltArousal +
            s(StandardizedTime, by=interaction(Group,Music)) + 
            s(StandardizedTime, Subject, bs="fs", m=1),
            AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M3d)
aicc_M3d <- AICc(M3d)
model_parameters(M3d)

--------------------------------------------------------------------------------
M4 = bam(Rating ~ Group*Music + Preference +
           s(StandardizedTime, by=interaction(Group,Music)) +
           s(StandardizedTime, Subject, bs="fs", m=1),
           AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M4)
aicc_M4 <- AICc(M4)
model_parameters(M4)

--------------------------------------------------------------------------------
M5 = bam(Rating ~ Group*Music + Preference +
           s(StandardizedTime, by=interaction(Group,Music)) +
           s(Subject, bs="re", m=1),
           AR.start=AR.start, rho=0.99, method="ML", data=dat)
summary(M5)
aicc_M5 <- AICc(M5)
model_parameters(M5)

--------------------------------------------------------------------------------
                     # Plot Figures for model diagnostics
--------------------------------------------------------------------------------
M5 = bam(Rating ~ Group*Music + Preference +
           s(StandardizedTime, by=interaction(Group,Music)) +
           s(Subject, bs="re", m=1),
           AR.start=AR.start, rho=0.99, method="REML", data=dat)
summary(M5)
k.check(M5)
model_parameters(M5)

--------------------------------------------------------------------------------
# Plot Figure S2 for model diagnostics
# Autocorrelation plot
tiff(filename="FigureS2a.tiff", width=6, height=4.5, units="in", res=900)
par(mgp=c(2.5, 1, 0), mar=c(5, 5, 2, 1))
check_resid(M5, split_pred='AR.start', select=3, ask=FALSE)
mtext("A", side=3, line=-2, at=0.068, outer=TRUE, cex=1.166, font=2)
dev.off()

# Normal Q-Q plot
tiff(filename="FigureS2b.tiff", width=6, height=4.5, units="in", res=900)
par(mgp=c(2.5, 1, 0), mar=c(5, 5, 2, 1))
check_resid(M5, select=1, ask=FALSE)
mtext("B", side=3, line=-2, at=0.068, outer=TRUE, cex=1.166, font=2)
dev.off()

# Combine Figures S2a and S2b
tiff("FigureS2.tiff", width=12.5, height=4.5, units="in", res=900)
FigA <- rasterGrob(as.raster(readTIFF("FigureS2a.tiff")), interpolate=TRUE)
FigB <- rasterGrob(as.raster(readTIFF("FigureS2b.tiff")), interpolate=TRUE)

grid.arrange(FigA, FigB, ncol=2, widths=c(1, 1), padding=unit(1, "lines"))
dev.off() 

--------------------------------------------------------------------------------
# Plot Figure 1 
# Panel A: Comparison of Musicians and Nonmusicians
tiff(filename="Figure1A.tiff", width=12, height=9, units="in", res=900)
par(mfrow=c(2, 2), mar=c(4, 4, 2, 1), oma=c(0, 0, 2, 0)) 

# Plot 1: Musicians vs. Nonmusicians in Tonal Music
plot_smooth(M5, view="StandardizedTime", cond=list(Group="Musicians", Music="Tonal music"), 
            col="#0072B2", lwd=3, lty=1, rug=FALSE, shade = TRUE, hide.label=TRUE,
            main=" ", xlab="Standardized Time (%)", ylab="Expectancy Rating", las=1,
            ylim=c(-100,100), cex.lab=1.4, cex.axis=1.3, font.lab=2)
plot_smooth(M5, view="StandardizedTime", cond=list(Group="Nonmusicians", Music="Tonal music"), 
            col="#009E73", lwd=3.5, lty=2, rug=FALSE, shade = TRUE, add=TRUE, xpd=TRUE)
legend('bottom',
       legend=c("Musicians Tonal music", "Nonmusicians Tonal music"),
       lty=c(1,2), lwd=c(3,3), col=(c("#0072B2","#009E73")), seg.len=3.6,
       bty='n', cex=1.4, ncol=1, xpd=TRUE)

# Plot 2: Musicians vs. Nonmusicians in Atonal Music
plot_smooth(M5, view="StandardizedTime", cond=list(Group="Musicians", Music="Atonal music"), 
            col="#D55E00", lwd=3, lty=1, rug=FALSE, shade = TRUE, hide.label=TRUE,
            main=" ", xlab="Standardized Time (%)", ylab="Expectancy Rating", las=1,
            ylim=c(-100,100), cex.lab=1.4, cex.axis=1.3, font.lab=2)
plot_smooth(M5, view="StandardizedTime", cond=list(Group="Nonmusicians", Music="Atonal music"), 
            col="#CC79A7", lwd=3.5, lty=2, rug=FALSE, shade = TRUE, add=TRUE, xpd=TRUE)
legend('bottom',
       legend=c("Musicians Atonal music", "Nonmusicians Atonal music"),
       lty=c(1,2), lwd=c(3,3), col=(c("#D55E00","#CC79A7")), seg.len=3.6,
       bty='n', cex=1.4, ncol=1, xpd=TRUE)

# Plot 3: Group Difference in Tonal Music
plot_diff(M5, view="StandardizedTime", hide.label=TRUE, las=1, rm.ranef=TRUE, 
          comp=list(Group=c("Musicians", "Nonmusicians")), 
          cond=list(Music="Tonal music"), main=" ",
          xlab="Standardized Time (%)", ylab="Difference in Rating",
          col="black", lwd=3, ylim=c(-100, 100), cex.lab=1.4, cex.axis=1.3, font.lab=2)

# Plot 4: Group Difference in Atonal Music
plot_diff(M5, view="StandardizedTime", hide.label=TRUE, las=1, rm.ranef=TRUE, 
          comp=list(Group=c("Musicians", "Nonmusicians")),
          cond=list(Music="Atonal music"), main=" ",
          xlab="Standardized Time (%)", ylab="Difference in Rating",
          col="black", lwd=3, ylim=c(-100, 100), cex.lab=1.4, cex.axis=1.3, font.lab=2)
mtext("A", side=3, line=-2, outer=TRUE, cex=1.5, font=2)
dev.off()

# Panel B: Comparison of Tonal and Atonal music 
tiff(filename="Figure1B.tiff", width=12, height=9, units="in", res=900)
par(mfrow=c(2, 2), mar=c(4, 4, 2, 1), oma=c(0, 0, 2, 0)) 

# Plot 1: Tonal vs. Atonal music in Musicians
plot_smooth(M5, view="StandardizedTime", cond=list(Group="Musicians", Music="Tonal music"), 
            col="#0072B2", lwd=3, lty=1, rug=FALSE, shade = TRUE, hide.label=TRUE,
            main=" ", xlab="Standardized Time (%)", ylab="Expectancy Rating", las=1,
            ylim=c(-100,100), cex.lab=1.4, cex.axis=1.3, font.lab=2)
plot_smooth(M5, view="StandardizedTime", cond=list(Group="Musicians", Music="Atonal music"), 
            col="#D55E00", lwd=3, lty=1, rug=FALSE, shade = TRUE, add=TRUE, xpd=TRUE)
legend('bottom',
       legend=c("Musicians Tonal music", "Musicians Atonal music"),
       lty=c(1,1), lwd=c(3,3), col=(c("#0072B2","#D55E00")), seg.len=4.0,
       bty='n', cex=1.4, ncol=1, xpd=TRUE)

# Plot 2: Tonal vs. Atonal music in Nonmusicians
plot_smooth(M5, view="StandardizedTime", cond=list(Group="Nonmusicians", Music="Tonal music"), 
            col="#009E73", lwd=3.5, lty=2, rug=FALSE, shade = TRUE, hide.label=TRUE,
            main=" ", xlab="Standardized Time (%)", ylab="Expectancy Rating", las=1,
            ylim=c(-100,100), cex.lab=1.4, cex.axis=1.3, font.lab=2)
plot_smooth(M5, view="StandardizedTime", cond=list(Group="Nonmusicians", Music="Atonal music"), 
            col="#CC79A7", lwd=3.5, lty=2, rug=FALSE, shade = TRUE, add=TRUE, xpd=TRUE)
legend('bottom',
       legend=c("Nonmusicians Tonal music", "Nonmusicians Atonal music"),
       lty=c(2,2), lwd=c(3.5,3.5), col=(c("#009E73","#CC79A7")), seg.len=4.0,
       bty='n', cex=1.4, ncol=1, xpd=TRUE)

# Plot 3: Music Difference in Musicians
plot_diff(M5, view="StandardizedTime", hide.label=TRUE, las=1, rm.ranef=TRUE, 
          comp=list(Music=c("Tonal music", "Atonal music")), 
          cond=list(Group="Musicians"), main=" ",
          xlab="Standardized Time (%)", ylab="Difference in Rating",
          col="black", lwd=3, ylim=c(-100, 100), cex.lab=1.4, cex.axis=1.3, font.lab=2)

# Plot 4: Music Difference in Nonmusicians
plot_diff(M5, view="StandardizedTime", hide.label=TRUE, las=1, rm.ranef=TRUE, 
          comp=list(Music=c("Tonal music", "Atonal music")),
          cond=list(Group="Nonmusicians"), main=" ",
          xlab="Standardized Time (%)", ylab="Difference in Rating",
          col="black", lwd=3, ylim=c(-100, 100), cex.lab=1.4, cex.axis=1.3, font.lab=2)
mtext("B", side=3, line=-2, outer=TRUE, cex=1.5, font=2)
dev.off()

# Combine Figures 1A and 1B
tiff("Figure1.tiff", width=25.5, height=9, units="in", res=900)
FigA <- rasterGrob(as.raster(readTIFF("Figure1A.tiff")), interpolate=TRUE)
FigB <- rasterGrob(as.raster(readTIFF("Figure1B.tiff")), interpolate=TRUE)

grid.arrange(FigA, FigB, ncol=2, widths=c(1, 1), padding=unit(1, "lines"))
dev.off()
--------------------------------------------------------------------------------
    # Repeated Measures Correlations between acoustic features and rating 
--------------------------------------------------------------------------------  
# Calculate rmcorr
data <- read_excel("data_4.xlsx")
names <- c("Subject","Group","Music")
data[,names] <- lapply(data[,names], factor)

# Split data into subsets based on Group and Music
musicians_tonal <-subset(data, Group == "Musicians" & Music == "Tonal music")
nonmusicians_tonal <-subset(data, Group=="Nonmusicians" & Music=="Tonal music")
musicians_atonal <-subset(data, Group=="Musicians" & Music=="Atonal music")
nonmusicians_atonal <-subset(data, Group=="Nonmusicians" & Music=="Atonal music")

# Helper function to compute repeated measures correlation
compute_rmcorr <- function(dataset, participant, measure1, measure2) {
  rmcorr(participant = get(participant), dataset = dataset,
         measure1 = get(measure1), measure2 = get(measure2),
         CI.level = 0.95, CIs = "bootstrap", nreps = 10000, bstrap.out = FALSE) }
--------------------------------------------------------------------------------
# Compute repeated measures correlations for each subset
##  Fluctuation Peaks
rmc_musicians_tonal1 <- compute_rmcorr(musicians_tonal, "Subject", "FluctuationPeaks", "Rating")
rmc_nonmusicians_tonal1 <- compute_rmcorr(nonmusicians_tonal, "Subject", "FluctuationPeaks", "Rating")
rmc_musicians_atonal1 <- compute_rmcorr(musicians_atonal, "Subject", "FluctuationPeaks", "Rating")
rmc_nonmusicians_atonal1 <- compute_rmcorr(nonmusicians_atonal, "Subject", "FluctuationPeaks", "Rating")
print(rmc_musicians_tonal1)
print(rmc_nonmusicians_tonal1)
print(rmc_musicians_atonal1)
print(rmc_nonmusicians_atonal1)

## RMS
rmc_musicians_tonal2 <- compute_rmcorr(musicians_tonal, "Subject", "RMS", "Rating")
rmc_nonmusicians_tonal2 <- compute_rmcorr(nonmusicians_tonal, "Subject", "RMS", "Rating")
rmc_musicians_atonal2 <- compute_rmcorr(musicians_atonal, "Subject", "RMS", "Rating")
rmc_nonmusicians_atonal2 <- compute_rmcorr(nonmusicians_atonal, "Subject", "RMS", "Rating")
print(rmc_musicians_tonal2)
print(rmc_nonmusicians_tonal2)
print(rmc_musicians_atonal2)
print(rmc_nonmusicians_atonal2)

## Key Clarity
rmc_musicians_tonal3 <- compute_rmcorr(musicians_tonal, "Subject", "KeyClarity", "Rating")
rmc_nonmusicians_tonal3 <- compute_rmcorr(nonmusicians_tonal, "Subject", "KeyClarity", "Rating")
rmc_musicians_atonal3 <- compute_rmcorr(musicians_atonal, "Subject", "KeyClarity", "Rating")
rmc_nonmusicians_atonal3 <- compute_rmcorr(nonmusicians_atonal, "Subject", "KeyClarity", "Rating")
print(rmc_musicians_tonal3)
print(rmc_nonmusicians_tonal3)
print(rmc_musicians_atonal3)
print(rmc_nonmusicians_atonal3)

## HCDF
rmc_musicians_tonal4 <- compute_rmcorr(musicians_tonal, "Subject", "HCDF", "Rating")
rmc_nonmusicians_tonal4 <- compute_rmcorr(nonmusicians_tonal, "Subject", "HCDF", "Rating")
rmc_musicians_atonal4 <- compute_rmcorr(musicians_atonal, "Subject", "HCDF", "Rating")
rmc_nonmusicians_atonal4 <- compute_rmcorr(nonmusicians_atonal, "Subject", "HCDF", "Rating")
print(rmc_musicians_tonal4)
print(rmc_nonmusicians_tonal4)
print(rmc_musicians_atonal4)
print(rmc_nonmusicians_atonal4)

## Novelty
rmc_musicians_tonal5 <- compute_rmcorr(musicians_tonal, "Subject", "Novelty", "Rating")
rmc_nonmusicians_tonal5 <- compute_rmcorr(nonmusicians_tonal, "Subject", "Novelty", "Rating")
rmc_musicians_atonal5 <- compute_rmcorr(musicians_atonal, "Subject", "Novelty", "Rating")
rmc_nonmusicians_atonal5 <- compute_rmcorr(nonmusicians_atonal, "Subject", "Novelty", "Rating")
print(rmc_musicians_tonal5)
print(rmc_nonmusicians_tonal5)
print(rmc_musicians_atonal5)
print(rmc_nonmusicians_atonal5)
--------------------------------------------------------------------------------
# Plot Figure 2
dat <- read_excel("C:/Users/Administrator/Desktop/data2.xlsx")
dat$Group <- as.factor(trimws(dat$Group))
dat$Music <- as.factor(trimws(dat$Music))
# print(levels(dat$Group))
dat$Subject <- factor(dat$Subject)
dat$Time    <- as.numeric(dat$Time)

# Standardize musical features
feat_names <- c("Fluctuation Peaks","RMS","Key Clarity","HCDF","Novelty","Rating")
dat[, feat_names] <- scale(dat[, feat_names])

# calculate mean rating
rating_summary <- dat %>%
group_by(Time, Group) %>%
summarise(mean_rating= mean(Rating, na.rm=TRUE), .groups="drop")

# plot 
plot_rating_acoustic <- function(acoustic_var, music_type, show_legend=FALSE,
                                 legend_position=c(0.535, 1.055)) { 
  dat_sub <- dat %>% filter(Music==music_type)
  rating_summary <- dat_sub %>%
    group_by(Time, Group) %>%
    summarise(mean_rating=mean(Rating, na.rm=TRUE), .groups="drop")
  
  acoustic_df <- dat_sub %>%
    group_by(Time) %>%
    summarise(acoustic=mean(.data[[acoustic_var]], na.rm=TRUE), .groups="drop")
  
  if (music_type == "Tonal music") {
    x_max <- 107
    x_breaks <- c(0, 107)
    col_mus  <- "#0072B2"
    col_non  <- "#009E73"
  } else {
    x_max    <- 146
    x_breaks <- c(0, 146)
    col_mus  <- "#D55E00"
    col_non  <- "#CC79A7"
  }
  
  p <- ggplot() + 
    geom_line(data= rating_summary %>% filter(Group=="Musicians"),
              aes(x=Time, y=mean_rating, colour="Musicians' Rating", 
                    linetype="Musicians' Rating"), linewidth=1) +
    geom_line(data=rating_summary %>% filter(Group=="Nonmusicians"),
              aes(x=Time, y=mean_rating, colour="Nonmusicians' Rating", 
                    linetype="Nonmusicians' Rating"), linewidth=1) +
    geom_line(data=acoustic_df,
              aes(x=Time, y=acoustic, colour="Acoustic Feature", 
                    linetype="Acoustic Feature"), linewidth=1) +
    scale_colour_manual(name="", 
      breaks=c("Musicians' Rating", "Nonmusicians' Rating", "Acoustic Feature"), 
      values=c("Musicians' Rating"=col_mus, "Nonmusicians' Rating"=col_non, 
               "Acoustic Feature"="black"), 
      labels=c("Musicians' Rating"="Musicians' Rating","Nonmusicians' Rating"=
               "Nonmusicians' Rating", "Acoustic Feature"=acoustic_var)) +
    scale_linetype_manual(name="", 
      breaks=c("Musicians' Rating", "Nonmusicians' Rating", "Acoustic Feature"), 
      values=c("Musicians' Rating"="solid", "Nonmusicians' Rating"="dashed", 
               "Acoustic Feature"="solid"),
      labels=c("Musicians' Rating"="Musicians' Rating", "Nonmusicians' Rating"=
               "Nonmusicians' Rating", "Acoustic Feature"=acoustic_var)) +
    
    labs(x="Time (s)",  y=expression(bolditalic(z)~bold(score)), title=NULL) +
    scale_x_continuous(limits=c(0, x_max), breaks=x_breaks) +
    scale_y_continuous(limits=c(-3, 5), breaks=c(-3, 5)) +
    theme_linedraw() +
    theme(axis.title.x=element_text(size=14, face="bold"),
          axis.title.y=element_text(size=14, face="bold"),
          axis.text.x=element_text(size=12),
          axis.text.y=element_text(size=12),
          panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(), 
          plot.margin=margin(t=12, r=6, b=12, l=6) )
  
  if (show_legend) {
    p <- p + theme(
      legend.position=legend_position,
      legend.justification=c(0.5, 1),
      legend.direction="vertical",
      legend.text=element_text(size=14),
      legend.key.width=unit(3.5, "lines"), 
      legend.background=element_rect(fill=alpha("white", 0), colour=NA),
      legend.box.background=element_blank(),
      legend.margin=margin(t=0, r=4, b=4, l=4)
    )
  } else {
    p <- p + theme(legend.position="none")
  }
  
  return(p)
}

p_tonal <- (
   plot_rating_acoustic("Fluctuation Peaks", "Tonal music", show_legend=TRUE) /
   plot_rating_acoustic("RMS",         "Tonal music", show_legend=TRUE) /
   plot_rating_acoustic("Key Clarity", "Tonal music", show_legend=TRUE) /
   plot_rating_acoustic("HCDF",        "Tonal music", show_legend=TRUE) /
   plot_rating_acoustic("Novelty",     "Tonal music", show_legend=TRUE) )

p_atonal <- (
   plot_rating_acoustic("Fluctuation Peaks", "Atonal music", show_legend=TRUE) /
   plot_rating_acoustic("RMS",         "Atonal music", show_legend=TRUE) /
   plot_rating_acoustic("Key Clarity", "Atonal music", show_legend=TRUE) /
   plot_rating_acoustic("HCDF",        "Atonal music", show_legend=TRUE) /
   plot_rating_acoustic("Novelty",     "Atonal music", show_legend=TRUE) )

col_label <- function(txt="A", size=13, face="bold", x=0.5, y=0.5) {
  ggplot() + xlim(0, 1) + ylim(0, 1) + theme_void() +
    annotate("text", x=x, y=y, label=txt, fontface=face, size=size*1.5/3.5) + 
    theme(plot.margin=margin(t=4, r=0, b=0, l=0) )  }
header_A <- col_label("A", size=13, face="bold", x=0.5, y=0.6)
header_B <- col_label("B", size=13, face="bold", x=0.5, y=0.6)

left_col  <- header_A/p_tonal  + plot_layout(heights=c(0.015, 1))
right_col <- header_B/p_atonal + plot_layout(heights=c(0.015, 1))

combined_plot_2col <- (left_col |plot_spacer() |right_col) + 
  plot_layout(widths=c(1, 0.05, 1))
ggsave("Figure2.tiff",combined_plot_2col, width=12, height=20, dpi=900)

################################################################################